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GREEN’S FUNCTION-BASED TIME STEPPING FOR THE 
KURAMOTO-SIVASHINSKY INITIAL-BOUNDARY VALUE 

PROBLEM* 

L. VAN VEEN? 


Abstract. Both theoretical and numerical studies of the Kuramoto-Sivashinsky equation have 
mostly considered periodic boundary conditions. In this setting, the Fourier decomposition of the 
solution is central to theoretical ideas, such as renormalization group arguments, as well as to nu¬ 
merical solution, allowing for the construction of accurate and efficient time-steppers using standard 
pseudo-spectral methods. In contrast, fixed boundary conditions induce boundary layers and neces¬ 
sitate the use of non-uniform grids, usually generated by orthogonal polynomials. On such bases, 
numerical differentiation is ill-conditioned and can potentially lead to a catastrophic blow-up of 
round-off error. In this paper, we use ideas recently explored by Viswanath (J. Comput. Phys. 251 
(2013), pp. 414-431) to completely eliminate numerical differentiation and linear solving from the 
time-stepping algorithm. We use the Green’s function-based method to investigate elements of the 
Kuramoto-Sivashinsky dynamics over a range of five decades of the viscosity. 
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1. Introduction. Boundary conditions form an important part of models of 
continuous, space and time dependent processes. Prescribing continuous variables 
like velocity or concentration, or their derivatives, at material boundaries often leads 
to boundary layers, i.e. regions with steep gradients of the dependent variables. Pe¬ 
riodic boundary conditions, in contrast, often lead to more spatially homogeneous 
solutions. Schemes for numerical simulation must take the peculiar structure of solu¬ 
tions, induced by boundary conditions, into account. 

Periodic boundary conditions arise naturally, and exactly, in some geometries, 
but are also often used when the actual domain is practically infinite or no more nat¬ 
ural choice is available. It is straightforward to list a number of reasons why periodic 
boundary conditions are pleasant to deal with. Firstly, the natural choice of a grid for 
spatial discretisation is a regular grid, and the natural choice of a basis to expand the 
solution in is one consisting of sines and cosines. The Fast Fourier Transform (FFT) 
provides an efficient way to switch between grid and spectral representations of a 
solution. Each basis function satisfies the boundary conditions and is an eigenfunc¬ 
tion of any spatial differential operator. Exploiting these facts, we can use standard 
pseudo-spectral methods to simulate the system. Examples of pseud-spectral methods 
for semilinear partial differential equations, such as the Kuramoto-Sivashinsky (KS) 
equation, can be found, for instance, in Trefethen’s text book |T^. The essential steps 
can be summarised as 

1. Discretize time using an implicit method for linear terms and an explicit 
method for nonlinear terms. This leads to a periodic Boundary Value Problem 
(BVP) to solve for each time step. 

2. Formulate the BVP in Fourier space. The linear differential operator is rep¬ 
resented by a diagonal matrix, so that the solution can be written explicitly 
in terms of the Fourier transform of the nonlinear terms. 
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3. Use the inverse FFT to find the dependent variables and their derivatives on 
a regular grid, evaluate the nonlinear terms there and use the FFT to find 
their contribution to the BVP. 

On a fixed number of n grid points in each spatial dimension, this approach yields 
a time step that requires O(nlnn) FLoating point OPperations (FLOPs) for each 
spatial dimension and has a bound on the spatial discretization error of the form 
exp(—cn), for come positive constant c, provided that the solution is smooth. 

With fixed boundary conditions, some complications arise. The Fourier bases 
and regular grids are no longer optimal, as they leads to spurious oscillations known 
as Gibb’s phenomenon. Instead, theory prescribes the use of orthogonal polynomials 
and clustered grids generated by their zeros or extrema. These basis functions are 
not eigenfunctions of all spatial differential operators and do not usually satisfy the 
boundary conditions. The most straightforward way to address these issues is to 
adjust the above scheme as follows: 

2a Formulate the BVP in spectral space. The linear differential operator is 
represented by a matrix that is structured (e.g. triangular) but not diagonal. 

2b Delete one linear equation for each boundary condition and replace it by the 
boundary condition on the polynomial basis. 

3 Use an appropriate (inverse) fast transform to evaluate the nonlinear terms 
on the clustered grid. 

4 Solve the resulting linear problem. 

In step 2b, only a small error is incurred for well-resolved solutions. Furthermore, for 
many polynomial bases, fast transforms from grid point values to and from expansion 
coefficients are know. For Chebyshev bases, for instance, the FFT can be used, as 
explained by Trefethen [16] . who presents several explicit examples of this approach. 

The most important issue, however, is the fact that the algorithm now requires a 
solver for a large linear system which is ill-conditioned for high-order differentiation 
and fine grids. For differentiation on Chebyshev bases the condition number grows as 
n 2p , where p is the order of differentiation. The KS equation is of fourth order and 
thus the condition number in a naive application of the Chebyshev spectral algorithm 
does not yield a useful bound on the accumulation of errors. 

These are various ways to improve the conditioning of the linear system. It 
appears that the most successful approach combines two elements: preconditioning the 
linear system by a spectral integration matrix and expanding the highest derivative of 
the unknown functions, rather than the unknowns themselves, in a polynomial series. 
Muite |12] compares several subtly different variations of this approach on grids up 
to 10,000 points. The preconditioned system still has a condition number that grows 
algebraically with the number of grid points and is of no use when estimating errors. 
A detailed examination of the linear systems arising in the various forms of spectral 
integration by Viswanath m lead to the conclusion that some can yield results near 
machine accuracy in spite of the bad conditioning. This good accuracy hinges on the 
cancellation of errors and the careful implementation of the boundary conditions. 

In the current paper, we follow a different, more radical approach, which was 
also investigated by Viswanath [18]. We solve the BVP analytically using Green’s 
function, thereby eliminating the need for numerical differentiation and linear solving 
altogether. Instead, we must construct a proper quadrature, that has an exponentially 
small error in spite of the finite differentiability of Green’s function. 

We show, that a combination of barycentric re-interpolation to sub grids and 
local Clenshaw-Curtis quadrature leads to an accurate time-stepping algorithm that 
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is stable to over resolution and requires only moderate computational resources to 
handle grids with tens of thousands of points. This enables us to generate numerical 
solutions to the KS Initial Boundary Value Problem (IBVP) for small values of the 
viscosity or, equivalently, on very large domains. The goal of those computations is 
to generate high quality data on the statistical behaviour of the solutions to the KS 
equation. In particular, we would like to find evidence for, or against, a conjecture 
put forward by Yakhot [22], which states that certain average properties of transient 
behaviour are described by the exponents of the Kardar-Parisi-Zhang equivalence class 
[§]. To the best of our knowledge, no other time-stepping algorithms with spectral 
accuracy have been developed for this purpose 

2. The KS IBVP. We consider the following IBVP, originating in the work of 
Kuramoto m and Sivashinsky [15] : 

Ut + UUx 'U'xx L / 'U j xxxx = 0 ( 2 . 1 ) 

u(-l) = l, u(l) = r, u xx (- 1) = u xx ( 1) = 0 (2.2) 

Other boundary conditions can of course be considered, fixing for instance u and u x 
at the boundaries. In that case, the time-stepping algorithm described here does not 
change, but Green’s function does, i.e. the analysis presented in Appendix [A] must be 
adjusted. In the literature, the KS equation is often presented in t he s caling u = yfvu , 
x = {x + \)/y/v and t = t/v, which yields an equation identical to ( |2.1[ ) but with v = 1, 
considered on [0,T] = [0, 2/y / IZ|. This scaling is used, for instance, in the numerical 
work of Wittenberg and Holmes [20] and the theoretical work of Galaktionov et a 1. 
0. The former work provides an overview of typical dynamics generated by the KS 
equation and a rich reference list. In the latter work, the existence of a bounded 
solution for any finite time of this IBVP is proven. 

Another common guise is the integral formulation 

ht = 2 (^v) hxx vhxxxx (2.3) 

where —h x = u. This form is often used in the physics literature when considering 
the KS equation as a model for interface growth, for instance in Refs. [2211]. For 
our time-stepping algorithm it is more convenient to solve for u on the fixed spatial 
domain and consider the viscosity v as control parameter. We then typically observe 
a transition from equilibrium at large values of v to time-periodic motion and finally 
spatio-temporal chaos as v is decreased. The spatio-temporal chaos exhibits a form 
of extensivity, as demonstrated in section [bj 

3. Time-stepping based on Green’s function. In this section, we will first 
describe the Green’s function-based algorithm on a high level, before specifying and 
justifying our choice of grids, interpolation and quadrature rules. 

It is convenient to consider deviations from a linear profile: 

Vt + vv x + R v + (j)v x + R (j> + v xx + v 

Vxxxx 0 (3.1) 

u(-l) = u(l) = i> xx (-l) = v xx (1) = 0 (3.2) 

T — l 

u = v + (f), 4>{x) = l + R{x + !),!?= ——— (3.3) 


which yields an IBVP with homogeneous Dirichlet boundary conditions. 
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A 

k = 1 k = 2 k = 3 k = 4 

0=1 

h 

ot k 1 

Pk 1 

0 = 2 

2/i/3 

a k 4/3 —1/3 

Pk 2 -1 

0 = 3 

6h/ 11 

a k 18/11 -9/11 2/11 

Pk 3 -3 1 

o = 4 

12/i/25 

a k 48/25 -36/25 16/25 -3/25 
Pk 4 -6 4 -1 


Table 3.1 

Table of coefficients of SBDF formula ( 3.J\ ) for orders up to o = 4. 


3.1. Reformulation as a linear BVP. The first step is to turn the IBVP into 
a linear BVP for every time step. This is done by the use of an implicit-explicit time 
discretization. In particular, we use a Semiimplicit Backward Differentiation Formula 
(SBDF) j2] to obtain 


Cv k+ i = a s v k+ i- s + A ^ Psf{v k+ i- s ) (3.4) 

S=1 S=1 

where A is a constant that depends on the time step size h and on o, the order of the 
SBDF formula. The sub scripts denote the approximate solution at different times. 
We have introduced the linear operator £ and the part of the BVP that is treated 
explicitly, /, according to 


£ = 1 + A R + Ad xx + Avd 

xxxx 

f(v) = -vv x - cf>(x)v x - Rcf>(x) 

The definitions of A and the coefficients a and /3 are listed in Table [XT 


(3.5) 

(3.6) 


The time-discretized equation (3.41, together with boundary conditions (3.2), con¬ 
stitutes a linear BVP since all quantities on the right-hand side are known explicitly. 
The solution can be written in terms of Green’s function as 


v k+ i = ^ a sG * v k+1 _ s + A ff s G * f(v k+ i_ s ) 


(3.7) 


S= 1 


where the star denotes the convolution 


S=1 


i 

(G * v)(x) = J G(x, y)v(y) d y 
y=i 


This expression is not suitable for numerical quadrature because the second convolu¬ 
tion contains derivatives of the unknown function. We apply integration by parts to 
obtain 


v k +i 


^(a s + /3 s AR)G*v k+ i- s + A ^ /3 S DG* 

S=1 S =1 


(\ v l+ 1- 8 + (jyuk+i-s') - ARG*</) 
= G * h + DG * I 2 + J (3.8) 
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where DG denotes the derivative of Green’s function with respect to the variable of 
integration of the convolution and, for later convenience, we have introduced short¬ 
hand notation for the terms that appear in the convolutions with G and DG and the 
constant term. Now the only derivative remaining can be computed analytically so the 
problem of the bad conditioning of numerical differentiation has been eliminated. In 
its place, we face two new challenges. Firstly, we must compute Green’s function and 
cast it in a form suitable for numerical evaluation. In Appendix [Aj appropriate ex¬ 
plicit expressions are derived. Secondly, we must accurately compute the convolutions 
in the knowledge that the integrand is only once or twice continuously differentiable. 
For this end, we will use polynomial re-interpolation followed by standard quadrature. 


3.2. Re-interpolation and quadrature. When using classical, global quadra¬ 
ture rules, the finite differentiability of the integrands would lead to a fixed rate 
of convergence, much like a standard finite difference method for solving the IBVP 
( 2.1|2.2 ) would have given us. We can solve this issue by using separate quadratures 
for the sub domains [— l,x] and [cc, 1] on which G(x,y) is smooth. This will require 
reinterpolation of the integrand onto suitable, non-uniform sub grids. 

Let us denote the global grid by aq, * = 0,..., n and the left and right sub grids for 
the i th global grid point by x h /\ j = 0,, n L ’* and x*’ 1 , j = 0,..., n R ’ 1 . In addition, 


we will denote the approximations of a quantity a(x) on these grids by a, a L,i and a R ’ 1 . 
Then the convolutions with G and DG are represented by matrix-vector products 


b = 

G * a —» 

b = Ma 



n L,i 

n R,i 


Mij = 

Zj ° p ^PQ-^qj + Zj 

/^lR,i /^R,2 

^pq n qj 


p,q=0 

p.q =o 


b = 

DG * a —> 

b = Na 



n L,z 

n R,z 


Nij = 

\ 1 s~iL,i pL,i | \ 1 

Zj ^pq^qj + Zj 

r)R,i 
U pq n qj 


p,q =o 

p,q=0 



(3.9) 


where 


• B h,l is a ( n L,l + 1) x (n + 1) matrix of interpolation from the global grid to 
the i th left sub grid, and likewise for B R ’ Z , 

• G h,l and G L,Z are the ( n h,z + 1) x (n L,z + 1) diagonal matrices with the values 
G(xi,y ) and d v G(xi,y) take on the i th left sub grid, i.e. 


[Xi.X R ’ 1 ) 


— A 

'- r pq u pq 


8G(x, y) 
8y 


G% = 6 pq G(a 

and likewise for G R,Z and G R,Z and 

C L,l is a row vector of quadrature coefficients such that 


( Xi,Xp *) 


(3.10) 


C h,i a h,i as 



dx 


and likewise for C R,Z . 

3.3. Explicit expressions. Here, we will give explicit expressions and assess 
the expected degree of accuracy for the combination of barycentric interpolation on 
closed Chebyshev grids and Clenshaw-Curtis quadrature. 
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We will take both the global grid and the sub grids to consist of Chebyshev points 
of the second kind, i.e. 

:Vi = cos ( — ), i = 0,..., n 
\n J 

x y = (?/> L ’ l ocos) J = 0,.. .,n M with t/j L ’ l (x) = -1 + *(£* + 1)0+ 1) 


- R,i = O r ’ 1 o cos) ( ), J = 0,..., n R ’ 1 with ^0) = ^(1 - Xi)(x - 1) + 1 


JTT 


1 


(3.11) 


An explicit expression for the interpolation matrix can then be found in Berrut and 
Trefethen [3|. It is given by 


r>L,i _ 


( y W fc QQ-a: g ) \ 

Vfc=o w q {x h y -x k )j 


with Wq 


1/2 for q = 0 

< (-1) 9 for 0 < q < n 
(—1)"/2 for q = n 


(3.12) 


and likewise for £? R ’®. The extremal points of the local grids coincide with global grid 
points. This leads to a singularity in the formula above, so that we must separately 
specify that 


B 0 q ~ V 


B 


Ij,1 


= 5, 


qn 


&Oq 


B 


R,1 

n R q 



The Clenshaw-Curtis row vector is a product of a row vector of quadrature weights 
with a matrix representing the discrete cosine transform: 


Cp’ ('J'i T 1) dp ^ | d k F k p 
k =o 


where a p 






1/2 ifp = 0 1 n L ’ i 
1 otherwise 


' 1 if k = 

< 2/(1 — 4fc 2 ) if k is 

0 if k is 

1 ( kpir \ 

—r cos -r 

n L ’* V n L ’ 1 J 


0 

even and k > 0 
odd 


(3.13) 


and likewise for C R,X . 

The number points in the local grids, n L ’ 1 and n R should be as least as large as 
the number of points in the part of the global grid they span, i.e. n h ’ 1 > n + 1 — i and 
n R ’ 1 > i + 1. Under this condition, the error of interpolation from the global to the 
local grids is negligible as compared the the error of interpolation of the solution onto 
the global grid. Importantly, barycentric interpolation is stable under over resolution. 
If we set nk’ 1 = n R ’* = n+ 1 in the implementation described below, the results remain 
the same, at least up to machine precision. 


4. Implementation. The task of time-stepping the IBVP can now be split up 
into two steps. In the first step, we compute the quadrature matrices M and N for 
given i/, A and n and a specific choice of the number of points in the local Chebyshev 
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grids. In the second step we iterate a SBDF by simply computing matrix-vector 
products. Both tasks can easily be performed in parallel by distributing the points on 
the global grid, and the corresponding rows of the quadrature matrices, over processes. 
Here, we describe an MPI-based implementation, which has the advantage that the 
memory requirements for storing the quadrature matrices can be reduced. In the 
following, we label the processes p = 0,..., P — 1, and process number p will compute 
the solution on n p points with indices i p through i p . We further assume that P « n so 
that the work load can be almost evenly balanced, i.e. n p as n/P for p = 0,..., P — 1. 

Algorithm [l] shows a pseudo-code for the construction of the quadrature matrices. 
The most costly steps in the main loop over global grid points are number 4, the com¬ 
putation of the Clenshaw-Curtis quadrature weights, and number 6, the computation 
of the interpolation matrix. Assuming that the total number of local grid points is of 
the same order as the number of global grid points, these steps have a FLOP count of 
0(n 2 ), bringing the total for this algorithm to 0(ro 3 ). In the parallel implementation 
each processor will thus handle 0(n 3 /P). 


One more remark about the quadrature algorithm must be made. If the order of 
the global grid is large, it may happen that it coincides with points on local grids up 
to machine accuracy. When using standard double precision arithmetic, for instance, 
this start to happen for global grid orders upward from 20,000. This leads to divisions 
by zero when evaluating (3.121. This can simply be circumvented by detecting overlap 
up to finite precision of the grids and replacing each row of B that holds the coefficients 
of interpolation onto a overlapping local grid point by a row of zero elements and a 
single element equal to unity, just like for the extremal points of the local grid. It 
is a remarkable fact, explained in detail by Higharn [7j, that the evaluation of the 
elements of B is otherwise stable, in spite of the small denominators. In our test of 
the accuracy of the numerical quadrature in section [5j we replace rows corresponding 
to local grid points that are closer to global grid points than 10 -13 and the resulting 
computation is stable up to a global grid order of at least 74,000. 


Algorithm [2] describes the time-stepping. The loop over time steps includes the 
computation of the integrands in equation ( |3.8| ), which takes O (n p ) FLOPS on proces¬ 
sor p , an MPI all-to-all communication of O(n) elements and a matrix vector product 
with the elementary FLOP count of O (n p n). If the MPI routine takes a similar 
amount of time to complete as O(n) FLOPS, it is obvious that the communication 
time will be negligible as long as n p » 1, as we have assumed. 


Figure [8] shows the wall time taken for building the quadrature matrices and 
time-stepping 4000 times for grid sizes 1000 and 2000. This test was run on a clus¬ 
ter computer, on a node with 24 AMD Opteron 2.2GHz processors, 32Gb of RAM 
memory, 512Kb of cache and a QDR InfiniBand connection. Clearly, Algorithm [I] 
scales linearly as predicted by the FLOP count since it does not involve any commu¬ 
nication. Time-stepping according to Algorithm [2j on the other hand, shows approx¬ 
imately linear speed-up up to 9 processors only. This results depends on architecture, 
on a machine with two quad-core CPUs, for instance, the speed-up saturates at 2 
processors. There are two possible reasons for this limitation. Firstly, the execution 
of the all-to-all communication depends on the details of the hardware and the MPI 
implementation and is highly sensitive to concurrent use by other processes of cache 
memory. Secondly, the all-to-all communication is blocking and the execution time 
for the partial matrix-vector product may vary over processes. 


Lastly, we turn to the nontrivial question of initializing the time-stepping. For 
SBDF orders two and up, we need the solution at previous time instants. We propose 
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ALGORITHM 1: Computation of the quadrature matrices and constant term 

Input: viscosity v, time-stepping parameter A, boundary values l and r and global 
and local Chebyshev grid orders n, n L ’*, n R,1 : i = 0, ..., n. 

Output: (n + 1) x (n + 1) quadrature matrices M and N and (n + 1) vector J, 
distributed by rows over P processes. 

1 for i = if ... if do ’/, loop over rows stored on processor p 

2 .. . ' ' 

Set x b l , j = 0, ..., n L ’ 1 and x*’ 1 , j = 0, ..., n R ’* according to (3.11 1 . 


Compute the Clenshaw-Curtis row vectors according to (3.131 and set 

L s~iL,i 

a <— G 


b L <- C L 


a <— C 


C K 


Multiply by Green’s function or its derivative, i.e. set 




p—o 


\ 1 r L syL ,i 
2_j °p^pq 
p—0 


■ £ 

p=0 


R j R 

CLp\j pg 0 


tR^R,i 

Zj °p Kj pq 

p— o 


Compute the interpolation matrices according to (3.12) and set 


Mij - £ a pBpj + £ Nij <- £ + £ & r R r ; for j = 0,..., 

p—0 p—0 p—0 p—0 


De-allocate al, Pj, o R , & R . 

Compute Ji = 

9 return A7 P , A p and J p , i.e. rows if through if of the quadrature matrices and the 


constant term. 


ALGORITHM 2: Green’s function based time-stepping 


Input: time step size h, number of inner and outer iterations N, N 0 , o initial points 
i>—o+i,... Vo, boundary values l and r. 

M p , N p and J p are stored on processor p. 

1 for j = 1,..., N 0 do 

2 

3 

4 

5 

6 


MPI scatter Vk-o+i —> v°_ 0+i ■ ■ ■ f° r * = 1, 

for k = 1,..., IVi do 


Compute the integrands If and If of (3.8 1 . 
MPI gather the integrands: lf }2 , • ■ ■, If 2 ~* 
Set 


7l 2- 


/ root to all 


•/. all to all 


P 

u k +1 


MPI gather vf .,... vf 1 
Root: output Vk + </>■ 


■ v k . 


% all to root 


four possible solutions: 

1. Using the first order SBDF formula with a small time step. This is a com¬ 
monly used method to seed backward differentiation formulae, and is em¬ 
ployed by Ascher et a 1. j^. For the Green’s function based time-stepper, 
however, this method has its limitation. If the time step is taken very small 
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Fig. 4.1. Left: wall time for the computation of the quadrature matrices according to Algorithm 
[7| (red) and for time stepping according to Algorithm^ (blue) using the first order SBDF. Shown is 
the wall time averaged over 20 trial for 4000 time steps with grid orders n = 1000 and n = 2000. 
The black line indicates linear scaling. Right: density plot of the solution computed in this test. The 
viscosity is set to v = 2 x 10 —4 , the step size to h = 10“ 5 and the initial condition is uq = 


for fixed viscosity, Green’s function approaches a delta distribution as its 
length scale, 6 _1 (see (A.10)) goes to zero. If this scale becomes comparable 
to the spacing of the local grid near its end points, of order n~ 2 , an instability 
can occur. 

2. Richardson extrapolation from lower order. For instance, we can seed the 
second order SBDF by approximating u(x, h) first by two first order steps of 
size h/ 2, giving u^ix, h), then by a single step, giving Uh{x,h) and finally 
setting u(x,h ) « 2uf l / 2 {x,h) — Uh(x,h). Similar expressions can readily be 
found for higher order seeding. This method has two disadvantages. Firstly, 
it relies on the cancellation of error terms , which is likely unstable for the 
higher order versions, so that round-off error is introduced. Secondly, for each 
step with a different value of A, we need to repeat Algorithm [I] which has 
order 0(n 3 ) complexity. Therefore, this method is only practical for small 
grids and low SBDF order. 

3. Using a known exact solution to the KS equation. An exact solution can be 
found, for instance, in Parkes and Duffy [13]. We use it in section[5]to test the 
accuracy of the SBDF formulae. Strictly speaking, this is not a solution to the 
IBVP, but for small enough viscosity the boundary conditions are satisfied 
far beyond machine accuracy. 

4. Growing a solution from a small perturbation to the zero solution. If we 
compose the perturbation out of eigenmodes of the linear part of (3.1), we 
can compute the solution backward in time under the assumption that the 
nonlinear term in negligible. The disadvantage of this method is that it 
requires a long time integration for the perturbations to grow to finite size. 
This strategy is used in our computation of the finite-size effects in section [6] 


5. Convergence and stability tests. We present two tests to evaluate the 
accuracy of the Green’s function based time stepping. The first test demonstrates the 
exponential convergence of the quadrature computed in each time step. We used the 
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Fig. 5.1. Error of the numerical approximation of the convolution with Green’s function. Test 
function [5.1\ ) is considered with 12, 120 and 1200 oscillations in [—1,1]. The solid line denotes the 
theoretically expected error of the interpolation and Clenshaw-Curtis quadrature. 


test function 


£(x, k) = -- . \ _ — -7T — 7: cos(27 xkx) — \ 

sv ’ ' 1 + sin(7r/cx) 2 2 v J 2 


(5.1) 


which satisfies the homogeneous Dirichlet boundary conditions (3.2) and describes 2k 
oscillations on the domain. A similar function was used by Trefethen to demonstrate 
the convergence of standard spectral methods for functions that can be analytically 
continued in a neighbourhood of [—1,1] in the complex plane |16j . The continuation 
of our test function has poles at ±iln(l + ^/2)/(kn) that determine an upper bound 
for the error of Lagrange interpolation on the global grid [3j, as follows: 


max |£(:r, k) — p n (x) | < C exp ( - 
x-e[—1,1] V 


a 

k U 


where C is a positive constant and p n the interpolant of order n. For k ^ 6, we find 
that a « 0.28. In fact, a = ln(l + V^)/^ up to corrections of order 0(1/A: 3 ). This 
interpolation error sets, in turn, an upper bound for the error in the Clenshaw-Curtis 
quadrature, the difference being a constant factor m 

Let £ represent the test function evaluated on the global grid, and let uj represent 
££ evaluated on the global grid. Then we measure the quadrature error 


e q — I £ A/c^loo 


as a function of the grid order n for fixed v, h and k. In the first test, we set 
iq = v = 2 x 10 4 , h = A = 10 -5 and k = 6 to generate a function somewhat similar 
to the final state of the simulation shown in figure fright). In the second test, we set 
V 2 = v = 2 x 10~ 6 , h = A = 10~ 7 and k = 60 and in the third rq = v = 2 x 10 —8 , 
h = A = 10~ 9 and k = 600. The rationale for this choice of parameter values is that 
we expect the typical spatial scale of variation of the solutions to decrease as y/v for 
small viscosity, as demonstrated in section [6j 

The three data sets collapse onto a single curve if we plot e q as a function of 
n/k. Along this curve, the quadrature error decreases approximately as exp(— an/k) 
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Fig. 5.2. Left: relative error in u after time-stepping over t/v = 8 for SBDF formulae of order 
up to four. The expected error is indicated by a solid line for each order. Right: initial and final 
solution in this test. 


as predicted. This results indicates that the main error introduced by the numerical 
evaluation of the convolutions is that of polynomial interpolation. We can make two 
further observations. Firstly, the method is stable to over resolution, as the error 
does not increase beyond the minimum around n/k = 120. Secondly, the minimal 
error appears to be determined by a build-up of round-off error in the matrix-vector 
product. Every increase in the number of grid points, and thereby vector elements, 
by a factor of ten leads to the same increase in the minimal error for n/k ^ 120. 

The second test illustrates the order of convergence of the SBDF formulae. For 
this end, we use the known exact solution to the KS equation mentioned above, which 
is given by 


'j(x, t) = c + 


15 VTl 
i9 7H; 


lltanh 3 ^-^=(x — ct — Xq)^ — 9tanh ^-^=(x — ct — xq) 


. (5 - 2) 

Here, c is the speed of a solitary wave connecting two constant solutions, Xq is its initial 
position and q = y/ll/76 is constant. Of course, w is is not an exact solution to the 
initial boundary value problem, but it approaches its left and right limit values at a 
rate of exp(—2 qd/y/v) a distance d away from the soliton. We have set v = 5 x 10 —5 , 
c = 1000 and Xq = —0.2 such that the variation of w at the boundaries, and the the 
magnitude of its second derivative there, are below 10 -25 for 0 ^ t/v ^ 8. A similar 
test was used by Anders [T] et a 1., but their viscosity is two orders of magnitude larger 
and, consequently, they were forced to consider time-dependent boundary conditions. 

Figure [5jleft) shows the error of time-stepping with SBDF formulae up to order 
four. In these tests, the Chebyshev grid order was fixed ton= 2000. For our choice of 
parameters, the exact solution has a complex singularity close to that of test function 
|5.1| with k = 6, and consequently we do not expect the interpolation error to play a 
role. The smallest error achieved is instead determined by the accumulation of round¬ 
off error, inversely proportional to the step size h , as indicated with a solid line. The 
initial and final condition in this test are shown in figure [5j right), which shows only 
the centre of the domain. 

Finally, we consider the stability of the time-stepping method. Figure[5]illustrated 
the stability for grid orders n = 1000 (left) and n = 2000 (right). The solution was 
initialized by a linear combination of eigen modes of the linear operator with small. 
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Fig. 5.3. Stability of the time-stepping algorithm for SBDF formulae of order up to four. The 
algorithm is stable below the solid lines connecting the symbols for each order and to the right of the 
dashed line. Left: grid order n = 1000, right n = 2000. The solid line at the top of the figures denotes 
the transition from imaginary to real-valued roots pO] of the eigenvalues of the linear operator. The 
dashed line indicates the point where y/v = 2 /n, i.e. the largest grid spacing is roughly equal to the 
typical length scale of the solution. See text for initial condition. 


randomly chosen amplitudes. A combination of viscosity and time step size was 
labeled stable if the integration proceeded up to t/v = 150 without blow-up. There are 
two boundaries to the stable regime. One lies close to the transition from imaginary 
to real roots of the eigenvalues of the linear operator given in |A.5[ If these roots 
are real-valued, Green’s function exhibits global oscillations, meaning that updates to 
the solution become dependent over an arbitrarily long distance in a single time step, 
which renders the step unstable. The other occurs when the typical length scale of 
the solution, expected to scale as *Jv, equals the maximal grid spacing, which is fixed 
in these experiments. 

6. Example computations: finite-size effects. To demonstrate the power of 
the time-stepping method based on Green’s function and Clenshaw-Curtis quadrature, 
we generated time series, seeded with random initial conditions of small amplitude, 
for a range of five orders of magnitude of the viscosity. Each time series extends up to 
t/v = 2000 and employed the SBDF formula of fourth order. After a transient time 
of about t/v = 150, the amplitude saturates and the dynamics is highly nonlinear. 
A fragment of each time series is shown in figure [6] For the smallest viscosity, v = 
1CT 7 , we have enlarged one tenth of the domain. The fact that the dynamics look 
qualitatively the same as that for v = 10 -5 is indicative of scaling behaviour, i.e. for 
small enough viscosity the solutions look similar in the scaled variables u , x and t 
introduced in section [2] 

Figure [6] shows the deviation from the time-mean solution near the left boundary 
in scaled variables to illustrate the dependence of the boundary layer thickness on 
viscosity. The fact that the curves approximately overlap indicates that the thick¬ 
ness of the boundary layer, in which the boundary conditions strongly influence the 
dynamics, is about \2^/v and is constant in scaled variables on the domain [0,T]. 

7. Discussion. We have demonstrated, that the Green’s function based method, 
in conjunction with the SBDF formulae, barycentric interpolation and Clenshaw- 
Curtis quadrature, is accurate, stable and reasonably fast for values of the viscosity 
as small as 10 -7 . To the best of our knowledge, no other time-stepping method for 
the KSIBVP has been tested for a viscosity this small or, equivalently, a domain this 
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large. 

Early studies of the KSIBVP, mostly in the physics literature, used finite-difference 
discretizations. Typical examples of such work are Manneville m and Sakaguchi M- 
Unfortunately, in this and similar work little description is usually given of the nu¬ 
merical methods, their accuracy and stability - Manneville’s work being a notable 
exception. There is no evidence that finite-difference-based results with a viscosity as 
low as 1CT 6 are reliable even for statistical analysis or curve-fitting exercise like that 
in the work cited above. 

In later work, various global and piecewise collocation methods were tested on 
the KSIBVP. Often, validation was only done for smooth, viscous solutions, like in 
Khater & Ternsah B|, who used spectral integration on a Chebyshev polynomial basis. 
Fornberg [5] applied a Chebyshev pseudospectral method, implementing the boundary 
conditions in real space, and computed chaotic solutions at v as 4 x 10 -4 . Piecewise 
collocation was tested for a viscosity of order O(10 -5 ) by Xu & Shu [23] and Anders 
et a 1. £T] . Based on careful testing and comparison to a priori error estimates, they 
could conclude that the spatio temporal chaos they observe numerically is a genuine 
property of the KS equation rather than a numerical artifact, but it remains unclear 
if their methods are suitable for time-stepping at smaller viscosity. 

The limiting factor of the Green’s function based method described here is the 
memory requirement, as it requires storing two (n + 1) x (n + 1) matrices - albeit 
distributes over processors - and the limited scaling of the parallelization. One pos¬ 
sible solution is to switch to piecewise Chebyshev grids to avoid excessive clustering 
of grid points near the boundaries. This will require a marginally more complicated 
procedure to compute the quadrature matrices and is work in progress. 

Acknowledgements. I would like to thank the Faculty of Engineering Science of 
Osaka University and the Japan Society for Promotion of Science for making possible 
the sabbatical visit during which much of this work was completed. 
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Appendix A. Computation of Green’s function. 

We are looking for Green’s function for the following linear BVP 

Cv = [1 + A R + Ad xx + Avd xxxx ] v = r t>(-l) = n(l) = v xx (-l) = v xx {l) = 0 

(A.l) 

where r will be set by all terms treated explicitly in the time discretization. With 
these boundary conditions, £, is symmetric, and has the following spectrum 

w% = sin(fc7ra;) A£ = 1 + A R — Air 2 k 2 + A^7r 4 /c 4 (A. 2) 

A| = 1 + A R — An 2 ^ + Att a u ^k — ^ (A.3) 

where the superscripts denote odd and even. We use the eigenfunction decomposition 
of Green’s function, given by 



GO, y) = Yj 


k =1 


mO )K(y) + w U x ) w Uy) 


K 


A? 


(A.4) 


The summations over the odd and even contributions proceeds in a similar fashion, 
so we will focus on the former. After factorising the denominator as 


Afc = Av(n 2 k 2 -p 2 + )(ir 2 k 2 -pi)-, p 2 ± = ± A 2 - 4Ai/(l + A R) 


(A.5) 


we can expand the summation as 
" /) _ 1 


z 

k= 1 




2Av(pl - p 2 _) ^ 


z 


COs(fc7T 0 — y\) COs(/C7T 0 — y]) 


7 r 2 k 2 — pi 


n 2 k 2 — pi 


-3? 

2 


2 A v{p 2 + -pi) 


cos(/c7r [a; + y]) cos(fc7r [a; + y]) 

7 T 2 k 2 — pi 7T 2 fc 2 — pi 

{/+ 0 — y) — f-( x ~ y) — f+ 0 + y) + f- 0 + v)}^) 


(A.6) 


where the ellipsis corresponds to some tedious manipulations of the sums to bring 
them into the form of the elementary inverse semi discrete Fourier transform 


f±(x) = 


1 


T 2 k 2 


■P± 


= z 


k =—00 


l 


r 2 k 2 — 


P± 


Airkx 


cos (p+ ~p±\x\) 
p± sin(p+) 


(A.7) 


As can be seen from this expr ession, Green’s function will exhibit global oscillations 

is positive so that at least one of p 2 + is real-valued. In 


A.5 


if the discriminant in Eq. 
that case, the resulting time-stepping scheme will be inaccurate and often unstable, 
as demonstrated in Sec. [5] We will therefore assume that 3(p+) # 0, which means 
that we impose an upper bound on the time step. 

Combining the odd and even contributions, we obtain Green’s function in the 
compact form 


G{x,y) = 


A i/Q(pi — pi) 


cos(2 p + -p+\x-y\) + cos{p+(x + y)) 
P+ sin(2p + ) 


(A.8) 
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In this form, Green’s function includes terms as large as exp(2Q(p + )) near the bound¬ 
ary, and p .|_, in turn, grows as l/y/v. This causes large cancellation errors near the 
boundaries. A more suitable form is 

G(x, y) = Q (e~ b \ x ~ v \ sin(2a — a\x — y\ — cj>) — e ~ 4b + b \ x -v\ s in(2a — a\x — y\ + (j>) 

—e~ 2b +b(x+y) s i n ( a ^ 2 ; + y) — (j>) + e ~ 2b ~ b ( x +y') sin(a(a; + y) + 4>)^j (A.9) 


where we have introduced the auxiliary parameters 

s= A 


4i/[l + A R] 

S-l/4 


(0 < S < 1) 


0-1/4 

a = 3?(p+) = 1 + Vs 

2 y/v 

0-1/4 /- 

b = %{p+) = 1 “ ^ 

h-Vs 

i + Vs 


9 = Arg(p + ) = arctan 

</> = Ar g (-^ 

\P+ sm(2 p+) 


(0 < 9 < 7t/4) 


= 2a — 9 — — + arctan 


sin(4a) 


e 4b — cos(4a) 
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1 — S ^/l — 2e -46 cos(4o) + e~ 8b 


(A.10) 


Where p + has been chosen to lie in the first quadrant in the complex plane. 

In this form, Green’s function no longer has exponentially large factors. It is 
immediately clear that, if the ratio between v and A is fixed, then the amplitude of 
G grows only as l/y'zl for small viscosity. However, there are still terms of 0(1) that 
cancel near the boundary. To avoid this, we rewrite Green’s function near x = 1 as 

G(x,y ) = Qsin(a(ir — 1)) ^cos(a(y + 1) + < p)(e b< ' x ~ y ~^ + e - fc G+y+ 2 ))_j 

cos (a(y + 1) - 4>)(e~ Kx ~ v) + e *>(*+?/- 2 ))J 
— 2Qsinh(6(a: — 1)) cos(a(x — 1)) [sin^y + 1) — (fi)e b( ' v ~ 1 V 

sin (a(y + 1) + (/)e _b ( y+3 )J (A.11) 


and use the latter form for numerical evaluation if 1 — x < 1/6. 
Green’s function satisfies 


(A.12) 


G(x, y) = G{y , x) G(-x, -y) = G(x, y) 

the latter property being a consequence of an S 2 symmetry of BVP (A.l I, namely 

( v , x, t) 


(- 


t) 


The BVP is equivariant under this reflection if, and only if, the original boundary 
conditions are, i.e. if r = — l in As a consequence of ( |A.12| ), we only have 

to evaluate Green’s function in on the domain 0 < x < 1, — x < y < x. Similar 
expressions are readily derived for G v (x,y) near the boundaries x = 1 and y = — 1. 
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Fig. 6.1. Fragment of the time series of u for five orders of magnitude of the viscosity. From 
top to bottom: v = 10“ 3 (n = 1000J, v = 10“ 4 (n = 3000 ), v = 10“ 5 (n = 10000 ), v 10“ 6 
(n = 24000J and v = 10 —7 (n = 60000J. Shown are the contours of u = y/vu for 900 ^ t/v ^ 1000. 
The initial condition in each case is a small random perturbation of the zero solution, the SBDF 
order is o = 4 and the time step is h/ 1 / = 5 x 10“ 4 . For v = 10 —7 , an enlargement of one tenth of 
the domain is shown so that the qualitative dynamics can be compared to that at v = 10 —5 . 
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Fig. 6.2. Root-mean-square of the departure of the solution from the mean profile near the 
boundary. The brackets denote averaging from t/v = 200 to t/u = 2000 and the scaled variables are 
u = yjvu and x = (x + 




